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Stochastic field equations represent a powerful tool to describe the thermal state of a trapped 
Bose gas. Often, such approaches are confronted with the old problem of an ultraviolet catastrophe, 
which demands a cutoff at high energies. In [l|] we introduce a quantum stochastic field equation, 
avoiding the cutoff problem through a fully quantum approach based on the Glauber-Sudarshan 
P-function. For a close link to actual experimental setups the theory is formulated for a fixed 
particle number and thus based on the canonical ensemble. In this work the derivation and the 
non-trivial numerical implementation of the equation is explained in detail. We present applications 
for finite Bose gases trapped in a variety of potentials and show results for ground state occupation 
numbers and their equilibrium fluctuations. Moreover, we investigate spatial coherence properties 
by studying correlation functions of various orders. 
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I. INTRODUCTION 

The enormous progress in experimental and theoreti- 
cal studies of ultracold quantum gases leads to a much 
deeper understanding of quantum many body physics 0- 
Hj]. In particular, Bose-Einstein condensation in traps 
enables us to investigate profoundly quantum statisti- 
cal phenomena in finite systems. In more recent experi- 
ments, spatial and temporal coherences of ultracold 
Bose gases are investigated in terms of correlation func- 
tions. Moreover, dynamics from nonequilibrium to equi- 
librium states [lfj and the spatial dependence of equi- 
librium density fluctuations [11| are considered. For a 
theoretical description of these phenomena it is desirable 
to consider nonzero temperature, finite size and - as no 
particle bath is present - a canonical description of the 
many-body quantum system. For most experiments, in- 
teraction between the atoms are of crucial importance. 
In systems with a Feshbach resonance, the interaction 
strength may even be tuned over a wide range - allowing 
to study ideal gases, too. In this article we deal almost 
exclusively with an ideal gas and comment on the inter- 
acting case briefly at the end of our paper. 

It is clear that many properties of ultracold gases are 
sensitive to temperature not least due to the different 
size of the condensed fraction of the gas. Moreover, for 
these finite systems with a fixed number of particles, 
it is preferable to base a theoretical description on the 
canonical ensemble rather than the usual grand canoni- 
cal ensemble. For an ideal gas, the best known example 
for the significance to choose the canonical ensemble is 
the fluctuation of the ground state occupation number 
(Nq-(No) 2 ) [l2j]- Clearly, it approaches zero as tempera- 
ture tends to zero in a canonical ensemble. Yet it is of the 
size of the average particle number N in a grand canon- 
ical description. A detailed study of the fluctuations of 
the ground state particle number in the microcanonical 
and the canonical ensembles arc given in [l3l - [l8j . Spatial 
correlations arc among other quantities of relevance for 



which a grand canonical formulation differs significantly 
from a canonical description [l9j . 

Recently, we presented a norm preserving stochastic 
field equation that describes the canonical state of an 
ideal Bose gas The resulting equation fully reflects 
quantum statistics and was given in [l[ in a representa- 
tion independent form. Therefore, it can be propagated 
in position space, meaning that it is not necessary to 
know cigenfunctions and -energies of the trapped atoms. 
While in [l[ we simply state the equation and display a 
few applications, we here want to elaborate on its deriva- 
tion, its numerical implementation and further applica- 
tions in much more detail. 

As we will show, spatial correlations and fluctuations 
can easily be calculated in position space. Due to the 
quantum framework, the equation - despite represent- 
ing a c-number field - docs not suffer from any cutoff 
problems which usually occur for classical field equations. 
Thus, numerical simulations can be performed with suf- 
ficient precision and reasonable effort. 

Let us relate our approach to previous stochastic equa- 
tions for the grand canonical ensemble. Note that most 
of these investigations are concerned with an interacting 
gas, while here, as a first step, we restrict ourselves to the 
case of the ideal gas. Since our theory can be formulated 
in position space, we clearly expect to be able to include 
interactions in a mean field sense in a next step. 

Detailed discussions of different approaches to (inter- 
acting) ultracold gases in terms of stochastic field equa- 
tions are given in [2(| [2l[ . Exact methods based on the 
positive P-rcprcscntation of the full density operator are 
used by Drummond and co-workers [22|. This theory 
brings about numerical challenges and is applied mainly 
to one-dimensional systems. More recently, also a 3D free 
gas at temperatures well below the transition tempera- 
ture has been discussed [23| . Other approaches [2(1 [23 - 
\2l\ can be seen as classical field methods, in which the 
lowly occupied energy levels must be treated in a differ- 
ent formalism or neglected, due to ultraviolet problems. 
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Davis and co-workers apply a projection on a subspace of 
highly occupied energy levels, using a cutoff [U], while 
Gardiner and co-workers separate the field operator in 
a highly and lowly occupied part [2(J [Hj]. We see our 
approach more in the spirit of Stoof and co-worker, who 
start from a path integral approach and arrive at a func- 
tional Fokker-Planck equation for the Wigner functional 
of the thermal field |26|. After additional simplifications, 
in their work the corresponding Langevin equation takes 
the form of a classical stochastic field equation as given 
by Hohenberg and Halperin [28J. 

Apart from these grand canonical approaches, the- 
ories for a fixed particle number have also been pre- 
sented more recently (29|-[32j . Xu and co-workers describe 
Bose gases with first order perturbation theory and in a 
Bogoliubov approximation [29|. An exact phase-space 
description for a finite number of particles is worked 
out by Korsch and co-workers. Exact evolution equa- 
tions for the Husimi-Q- and the Glauber-Sudarshan-P- 
distribution function based on a Bosc-Hubbard Hamilto- 
nian are derived [3Cj H3 ■ Stochastic methods in phase- 
space for the canonical ensemble of qubit and spin models 
have been presented by Drummond and co-workers [32j . 
While all these theories based on a canonical ensemble 
focus on different systems, e.g. lattices or spin models, 
we here aim to describe the full canonical thermal state 
of a Bose gas in any given trap. 

Our novel equation presented in [l[ has been found 
on the basis of the Glauber-Sudarshan-P-representation 
- being exact for an ideal gas. In this paper we aim to 
explain in more detail the derivation of the equation, its 
properties, its numerical implementation, and we show 
its power by discussing further applications. The pa- 
per is structured as follows: In the second section we 
show the theoretical background and the derivation of 
the stochastic field equation. The third section is de- 
voted to the discussion of some crucial properties of the 
equation. In the fourth section we demonstrate how we 
propagate the equation in energy and in position space 
numerically. Results of these numerical applications are 
shown in the fifth section, before we draw conclusions in 
the final section. 



II. STOCHASTIC FIELD EQUATION FOR THE 
CANONICAL ENSEMBLE 



In order to derive a quantum stochastic field equation 
for the canonical ensemble of an ideal Bose gas, our first 
aim is to express quantum expectation values for a fixed 
particle number in terms of c-number integrals. The de- 
sired stochastic field equation should enable us to calcu- 
late these expectation values. We use the corresponding 
Fokker-Planck equation to demonstrate the equivalence 
of stochastic and canonical quantum ensemble mean. 

We begin with the canonical density operator for N 



particles 



Pn 



Z 



Pn- 



N 



(1) 



As usual, the Hamiltonian H = ^ekd k dk can be ex- 

fc 

pressed in terms of occupation numbers of single-particle 
states (efc denotes the energy of single-particle state k). 
Furthermore, the canonical partition function is denoted 
by Zjv, with the usual = pf incorporating Boltz- 
mann's constant k and temperature T. Crucially, a pro- 
jector Pn = \{ n k})({nk}\ on the iV-particle sub- 

J2n k =N 

space appears, with the usual number states = 
1 n("l)™' t The operator P/v projects the expo- 

Un k \ k 

nential e~@ H on a subspace of N particles. In order to 
arrive at a c-number representation of quantum expec- 
tation values, one can express e~^ H from ([1]) as a mix- 
ture of coherent states u sing the (Glauber-Sudarshan) 
P-function representation |33( 

d^{z}P({z})\{z}){{z}\. (2) 

Here, we abbreviate the normalization factor with n = 

rj(e^ £fe — l) -1 , and products of coherent states with 
fc 

|{z}) = |zo)|^i) • • • \ z n) • • • ■ The appropriate measure is 

d/i{z} = n ( d{Rc " fc} ^ {lmzfc} S 

Of particular interest are correlation functions of quan- 
tum field operators. Using the P-representation 
we can determine canonical quantum expectation values 

with the help of ({z}\P N \{z}) = \z k \ 2 ) N e~ ? M ' 

fc 

(see [3J|) and obtain a weight function Wn({z}) = 



A 



-Ekf 



P({z}). With the latter, 



quantum expectation values can be written as integrals 
over c-numbers, 

(ajaj)iv = tr (a\a 3 p N ^j = — J dp{z}z*z ] W N - 1 {{z}) 



(3) 

with C = J d^i{z}W N ({z}). 

In this work we focus on spatial correlation func- 
tions. Therefore, we are interested in ensemble aver- 
ages for the bosonic field operator ip(x) = J2( x \ e k)dk 

k 

with \ek) the single-particle eigenstate corresponding to 
energy e^. First, the weight function is rewritten as 

WMM) = ^(fdxMxWfei-f^^Vpir,^) 
using if)(x) = z k(x\ek) which is the eigenvalue of the 

k 
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field operator when applied to the coherent state \{z}), 
i.e. ip(x)\{z}) = ip(x)\{z}). Further, we obtain 



(4) 



with the normalisation constant C = J d[i{ip}WN ({V'}) • 
It should be mentioned here that for quantum ex- 
pectation values of second (or higher) order one 
has to use weight functions of different N. For 
(^(x)^(x')ilj(x)if>(x')), for instance, the first order 
Wn-i({4>}) has to be replaced by Wn-2{{iP}) ~ f° r 
higher orders the result changes accordingly (k-th order: 
chose W N - k {{ip})). 

In order to calculate these expectation values we were 
able to construct a new, norm-preserving stochastic field 
equation [l[ . The equivalence of quantum statistical and 
stochastic ensemble mean is based on the fact that the 
weight functions W n (n = N,N — 1, . . .) turn out to be 
stationary solutions of the corresponding Fokkcr-Planck 
equation. 

In representation independent form NAx, t) = 
(x\ip{t))) and using Stratonovich calculus |35( it reads 
(h = 1 throughout) 



(S) dty) 



(A + i)H 



2 A VK\d£) 



\ip)dt 



KW) ■ (5) 



Here, the single particle Hamiltonian H = 2— + V(x) = 

X e fel e fe)( £ fc| appears, and a damping parameter A. Cru- 

k 

cially, temperature enters through an operator 



K = 



H 



(6) 



The noise increment \d£) is uncorrelatcd in space and 
time (x\d£(t)){d£(t)\x') = S(x - x')dt. Equation © is 
a norm preserving, non-linear stochastic equation that 
enables us to obtain canonical expectation values on av- 
erage. 

It is important to stress that the dependence of the 
stochastic fields \ip(t)) on the particle number N is in- 
directly given through their norm Af = which is 
preserved during propagation and thus determined by 
the initial condition. One might expect Af «a N to be 
a reasonable choice; however, matters are more delicate 
and a distribution of norms Af is required in order to ob- 
tain exact quantum expectation values. We devote the 
whole next Section IIIII to the relation between norm Af 
and particle number N. 

Sometimes it useful to express ([5]) as an Ito-stochastic 



equation, reading 

(i) d\i>) = - ( (a + i)H - a^S{k \ \m 




tr(AT) - 



K\i>)dt 



«-^4 p> 

The novel stochastic field equation can be easily solved 
numerically and used in different representations. An 
implementation in the position representation is par- 
ticularly useful as it can be applied to arbitrary ex- 
ternal potentials V(x). Based on the assumption of 
sufficient crgodicity, in applications the ensemble aver- 
ages are replaced by a long time limit (i/j1(x)iP(x'))n — 
t 

lim j [ ds {ib*(x, s)ip(x', s)} over a single realization 

ip(x, s) of the stochastic field equation. 

The term (A + i) consists of a phcnomcnological damp- 
ing term A, while "i" describes "real" dynamics. So, in 
a phcnomcnological manner, our equation can also be 
used to mimic the transition from a non-equilibrium to 
an equilibrium state. 

It is important to point out the crucial role of the op- 
erator K. Its meaning becomes most apparent when we 
omit the nonlinear terms in our equation ([5]) and think 
of H to represent H — fj, which also affects the operator 
K. One obtains an equation which gives us the thermal 
state of the grand canonical ensemble [HI 



d\ip) = -(A + i)H\ip)dt + V2AK\d£). 



(8) 



with \x the chemical potential used to fix the average 
particle number N. It should be mentioned here that the 
canonical equation ([5]) is not merely a normalized version 
of the grand canonical equation. The linear equation (JSJ) 
is closely related to the classical field equation for the 
thermal state of the grand canonical ensemble for any 
Hamiltonian energy functional ip*) [HI 



dip{x, t) 



IK ^ 



y/2AkTd£(x,t), (9) 



which has well-known ultraviolet problems. If one sets 
Jjp- = Hif, there is only one difference between classical 
© and quantum equation (J5}: the classical temperature 
kT is substituted by the operator K. As discussed in [2J], 
for instance, the classical equation satisfies the equiparti- 
tion theorem and therefore the infinite degrees of freedom 
of a field i/j(x, t) lead to an ultraviolet catastrophe. This 
can also be seen in equ. ([9]) looking at a formulation 
in position space: the white noise term which is uncor- 
relatcd in position induces arbitrarily high momentum 
kicks. However, if we consider the operator K = et >u _ x 
acting on the white noise \d£) in position space represen- 
tation, the former white noise becomes spatially corre- 
lated, and we no longer get unphysical momentum. In 
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{x\K\x') 



kT = 1.2hu 




{x\K\x') 



kT = 4.5ftw 




FIG. 1. Correlation function {x\K\x'} of the noise \d() — \JK of our stochastic field equation ([5j for a harmonic potential. 
We show the correlation function for temperature kT = 1.2ftu (left) and for a higher temperature kT = 4.5ftu (right). 



other words, proper quantum statistics requires the re- 
placement of spatially uncorrelated noise VkT\d£) in (J9]) 
by spatially correlated effective noise |g?£) = \J~K |<i£), as 
in ©. The correlation function of the effective noise 
(x\d((t)) (d((t)\x') = (x\K\x') dt is now a smooth func- 
tion in position space and shown in Fig. [T]for two differ- 
ent temperatures. 

In the low energy limit H <C kT, indeed, the operator 
tends to the classical temperature K w kT. Yet in the 
case of high energies H ^> kT, the operator tends to zero 
K — > 0. We see that in the quantum equations ([51 [5]) the 
required cutoff is built in. 



III. STOCHASTIC AVERAGE AND QUANTUM 
EXPECTATION VALUE 

In this chapter we want to investigate the relationship 
between quantum expectation values of arbitrary order 



djdk.-.di) ([3]) and those obtained with the stochastic 
field equation (SFE) ((z* ...z* z k ...zi)) SFE . For the follow- 
ing we introduce the norm J\f = E \ z k\ 2 = °f the 

fc 

wavefunction. The stochastic equation keeps the norm 
fixed, while the integral of the quantum expectation val- 
ues ([3|) extends over a range of values of the quantity 

E \ z k\ 2 - The norm stays constant and for any stochas- 

fc 

tic method ((1))sfe = 1- Hence, we can write the first 
order expectation values we obtain from the stochastic 
equation as 



(« 2: j>>sfe(A/', e ) 

/ d»{z}z* Zj 5(Af -EN 2 )WV-i(M) 

_ fc 

$dii{z}5(N~Y.WV)W N - 1 {{z}) ■ 

k 



(10) 



We find that this expression depends on the ground 
state energy eo- It can be seen easily by a sub- 



stitution z — > ze^ E(| / 2 that the simple identity 
e' 3e »((z*z J )) SF E(AAe-^,eo) = ((^»Sfe(JV, 0) holds. 
The relation between the average values given by the 
stochastic field equation and quantum expectation val- 
ues is now obvious 



(a\a j ) N = I dNP(N)((z*z 3 )) SFE (N,e )- 



with the norm distribution 



P{M) 



J dn{z}W N {{z}) 



(11) 



(12) 



in the integral (jH)) . 



Note that / dNP{N) 

J dn{z}W N ^({z})/ J dfi{z}W N ({z}) + 1. There- 
fore, and in order to get a better understanding of 
cqu. (|11|) . we have to investigate P{M) more closely. A 
lengthy calculation is given in appendix [XJ which results 



P(A0 



1 fc 

(N-l)\ ^cie-^'t^ 1 ) ' 

k 



(13) 



Well below the critical 



with c fc = II [ eH 

1=0 
l^k 

temperature it is sufficient to consider only the first term 



of the sum E c fc e 
fc 



because the remaining terms are 



smaller by a factor k ~ eP °W . Numerical investiga- 

tions of the factors c k showed that this approximation is 
justified in the temperature region considered here. Then 



5 



our distribution P(Af) is normalized to 



dAfP(N) 



k 



(14) 



Within this approximation we have a Poisson-likc distri- 
bution 



P(A0 ~ Af iN -V expf-e^AO- 
This distribution is centered around 

oo 

/ dMP{N)N 

W = ^5 = 

J dNP{N) 





(15) 



(16) 



and the standard deviation equals AAA = ^/Ne~^ e °. 
For large N the ratio of variance and mean goes to 
zero — > 0. Hence, in general it is sufficient 
to propagate equ. © with a single norm, i.e. we 
replace P(N) e pe °S(Af - 2Ve _/3eo ). As a conse- 
quence, from (|11[) we get the simple relation (ajeij)jv = 
e@ e ° ((z* Zj)}sFE(Ne~P e ° , cq) between canonical quantum 
correlation function and ensemble average of the stochas- 
tic field equation. As further elaborated upon in Section 
IVl it is for very high precision only, that it is necessary 
to take into account the full norm distribution. Let us 
also remark that for higher order expectation values these 
considerations apply similarly. 

It should be noted here that if the norm of our numer- 
ical simulation of the stochastic field equation Af should 
correspond to the particle number N, we have to choose 
£o = 0. However, it is also possible to choose a norm dif- 
ferent from the particle number according to (fTB"|) . This 
freedom is useful for an easier numerical implementation, 
as shown in Section [TV] It should also be mentioned that 
for temperatures above and near the critical tempera- 
ture, the approximations leading to expression (fl4j) are 
no longer valid and the exact norm distribution needs to 
be known. Still, using (| 14|) . for many quantities we ob- 
serve satisfying results even in this temperature region. 



IV. NUMERICAL IMPLEMENTATION 

In this section it is shown how one can solve the novel 
stochastic field equation numerically. The implementa- 
tions can be done in different representations, whichever 
is convenient for a given trap potential. If the single par- 
ticle eigencnergies are known (recall that we neglect inter- 
atomic interactions here) a propagation in energy space 
is the easiest way to solve our equation. In a box po- 
tential this is analogous to a propagation in momentum 
space. For a trap potential V{x) with unknown energy 
spectrum, one can solve the equation in position space 
directly, using a split operator method based on Fast 
Fourier Transformation. This implementation is more 



challenging but the program can be easily adjusted to 
any trap just by changing the single line of code where 
the potential is defined. 

If the eigenvalues for the given trap potential are 
known, a numerical implementation in energy space does 
not contain any challenging aspects. As already ex- 
plained above, it is more convenient to pass to position 
space, if we want to solve problems with arbitrary exter- 
nal potential, for which the eigenenergies are not known. 
We turn to a discretized formulation of equation ([7]) and 
use a position grid {xi} with some multi-index i. The 
part of the equation containing the single-particle Hamil- 

2 

tonian H = J~ + V(x) only can be propagated easily us- 
ing standard Fast-Fourier transformation methods [37| . 

The challenging aspect is the evaluation of matrix el- 
ements of the operator ([6]) appearing in (O, i.e. the de- 



termination of (xi\K\xj 



-+V(x) 



- + V(x) 



For an efficient implementation a Wigner-Weyl- 
approximation is helpful. The Wigner-Weyl correspon- 
dence |33| of the operator K in 2D-dimensional phase 
space is defined as 

K(x,p) = Jd D se^ s (x + %\K\x-i), (17) 

with the help of which we can obtain the matrix elements 
as 



{xi\K\xj) = J d D peW«-**) K (i 



■P 



(18) 

In order to simplify, we consider only zeroth order in h, 
and find the classical expression 



K(x,p) 



(V(a 



(19) 



This approximation will be justified by comparing the 
approximated numerical result with exact numerical cal- 
culations. 

Expanding the denominator in a geometric series we 
get 

(xi\K\ Xj ) « 77~t7t / dPpe***-** ( V ( Xl+X] 



(2tt)' 



£' 



-PnlVt 



2 I > 2 



(20) 



The momentum-integration can now be carried out 
without any difficulty which leads to the following ex- 
pression: 



(xi\K\xj) 



1 -0nV(^ 



V 



y^. 3_>_ — 

2f3n 2(/3n) 2 



(21) 
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FIG. 2. Correlation function {x\K\x') of the correlated noise of our stochastic field equation (J3J for a harmonic potential 
at temperature kT — 3hio. We show the correlation function obtained from the noise in Wigner-Weyl approximation (left) 
compared to the exact result (right). 




FIG. 3. Wigner function of a Bose gas of 100 particles trapped in a harmonic potential (left) and in the quartic potential 
V(x,y, z) — x 4 + y 4 + z 4 (right) above (top), at (middle), and below (bottom) the critical temperature. On the left hand side 
of both pictures we display a single realization with our stochastic field equation ([5]), on the right hand side a long-time average 
over 15000 time steps. 

Top: almost a Boltzmann distribution for T > T C (N); middle: a peak develops for T ~ T C (N); bottom: almost all particles 
condense in the ground state, the Wigner function reflects the ground state Wigner function T < T C (N). 



It is important to note that the convergence of the geo- 
metric series in (|20|) very much depends on the minimum 
of the potential energy. Indeed, rapid convergence can be 
assured by adding some (physically irrelevant) constant 
to V(x). However, as this shift also affects the ground 
state energy e , it also has a significant influence on the 
norm distribution of the wave function as can be wit- 
nessed in (fl6|) . Therefore, in practice one tries to identify 



some optimal shift that assures both, rapid convergence 
and a reasonable norm of the wave function. 

An efficient generation of (xi\K\xj) is given in equ. 
([21~j) . Now, in order to propagate equ. ([7]), it is neces- 
sary to think about the numerical implementation of the 
correlated noise \dQ = \/~K\d£), too. Obviously, the (dis- 
cretized) noise dQ = (xi\dQ must be generated in a way 
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that 



dQd(* 



(xi\K\xj)dt 



holds. In order to find an efficient way to obtain the 
noise, we start with equ. (|20|l . Based on the substitution 
p —> \/nj3p we get 



(xi\K\xA 



1 



D 



E 1 

n=l 
X 



dPpe 



1 i ad P ' i 
— — I a p—e 

n(3j P 2 



(22) 



For simplicity, the following derivation of an efficient 
noise generating method will be restricted to the one- 
dimensional case D = 1. For higher dimensions, the 
strategy follows the very same lines. 

.p(l£ i -Xj) ^2 

The first integral J dpe e~~ can be obtained 

from a Monte-Carlo integration of ((e l )) over a 

normally distributed real random variable r\\. For the 
second momentum integral in (|2"2"j) , note that the integral 

J dp^-e 1 e~~ can be seen as a radial integration 

in spherical coordinates according to 



/ 



p2 pC 

dp-e* 



dpp 2 cos 



p{xi 



x i)\ -si 
1 e 2 



4tt 



d p cos 



fnfi 



(23) 

Obviously, now the integral in equ. (|23p can be evaluated 
by a Monte-Carlo integration using three additional real 
Gaussian random numbers 772,773, and 774. We arrive at 
the following representation of the correlation function 
as an average over four normally distributed real random 
variables 771, . . . , 774: 



To proceed further we make use of the observation 
that the correlation function is very much concentrated 
along the diagonal x — x' (see Figs. I1I2[) . Based on 

this fact, we approximate V \ ~ y/V(xi)V(xj) 

and V ( X, + X] ] pa v(x l )+v(x J ) ^ reS p ec ^j V ely, in order to 



2 J 2 
achieve a factorization of the x-dependence on the right- 
hand-side in equ. ([2~4"[) . The final step is the introduc- 
tion of a family of independent complex Wiener incre- 
ments dKi n (with i = 1,2,3 and n = 1,2,3,...), i.e. 
dKi n dn* k = SijS n kdt. With their help we achieve an 
explicit expression for the correlated noise in terms of 
independent Gaussian (77) and Wiener (dn) random vari- 
ables, 



— ; (znpir W \ 




cine 



x [ dn 2n cos ( — ^=\/ V2+V3+ vl 



+dK 3n sin [ — ^=\/ 77I + r/l + r)\ 



(25) 



y/pj +pl+ pj{ Xi - Xj ) \ c _?i±4±A_ { Xi \K\ Xj ) 



The good agreement of the correlation function of the 
noise in Wigner-Weyl approximation from (|25j) and the 
exact correlation function can be seen in Fig. (5[ where we 
show on the left an average over numerical realisations 
from (|25p and on the right the exact correlation function 
for an isotropic harmonic oscillator at kT = 

To summarize this section, we are able to obtain the 
non-trivial noise and matrix elements of K with reason- 
able effort. The other terms of the equation can be prop- 
agated using Fast-Fourier transformation. Thus, we are 
able to propagate the whole stochastic field equation in 
D spatial dimensions. Results of such numerical simula- 
tions are presented in the following Section IVl 



(xAK\x~ 



1 



(27m0)* ^ 



w / x ■ 4- x ■ \ 



V 




(24) 



V. RESULTS OF NUMERICAL SIMULATIONS 

Our equation is here applied to a 3D Bose gas of fixed 
particle number trapped in various 3D potentials and dif- 
ferent quantities are calculated. The results are obtained 
in energy or position representation. The applicability of 
our equation is first shown by determining the Wigncr 
function W(x,p) = ± J dy e l ™ lp (x + \)${x - §)) of 
a Bosc gas. The propagation itself is performed in 3D 
position space, the figures show the Wigner function in- 
tegrated over the remaining four phase space coordinates. 
As already explained in Section [TV] we use a Wigner-Weyl 
approximation for the simulation of ([7]). 

In Fig. [3] we show the Wigner function for an ideal 
gas of 100 particles for different temperatures and for 
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an isotropic 3D harmonic potential (left) and the poten- 
tial V(x, y, z) — x 4 + y 4 + z 4 (right). The left columns 
in both figures show a snapshot of a single realization, 
while the right columns display an averaged result over 
15000 time steps. Clearly, to arrive at the equilibrium 
distribution, an average over many runs is required. The 
temperatures are chosen above, at, and below the crit- 
ical temperature T C (N) (for the finite N particle sys- 
tem). For the harmonic potential (left) one sees the ex- 
pected qualitative properties; when applied to the po- 
tential V(x, y, z) = x 4 + y 4 + z 4 (right) we see that our 
implementation can be easily used for arbitrary trap po- 
tentials with an unknown energy spectrum. 




+ SFE V(x,y,z) 


= x i + y i + z i 


x SFE V(x,y, z) 


= x 6 + y 6 + z e 


* SFE V(x, y, z) 


= x a + y 8 + z 8 


— th. 1. V(x,y,z) 


= x i + y i + z i 


— th. 1. V(x,y,z) 


= x e + y e + z 6 


■-■ th. 1. V(x,y,z) 


= x 8 + y 8 + z 8 



FIG. 4. Ground state occupation as a function of temperature 
from a simulation of an ideal 3D Bose gas of 200 particles 
trapped in the potentials V(x, y, z) — x 4 + y 4 + z 4 (plus signs), 
V(x,y,z) = x e + y 6 + z 6 (crosses) and V(x,y,z) = x 8 + y s + z 8 
(stars). We compare with the thermodynamic limit (full: x 4 , 
dashed: x 6 , dashed-dotted:a: 8 ). 

We next show that we also obtain good results quanti- 
tatively. First, we focus on the ground state occupation 
of the ideal gas. In [l| we showed the good agreement 
of our values with other canonical results for a box and 
a harmonic potential. In Fig. [4] results for the con- 
densed fraction in potentials with higher powers than 
the harmonic case are shown. We use the 3D poten- 
tials V(x, y, z) = x 4 -\- y 4 -\- z 4 (plus signs), V(x, y, z) ~ 
x 6 +y 6 + z 6 (crosses) and V(x, y, z) = x s + y s + z s (stars) 
also to illustrate the flexibility of the numerical imple- 
mentation of our equation. We compare with the results 
obtained in the thermodynamic limit (full line, dashed, 
dashed dotted); clearly, finite size effects are visible. 

Spatial correlation functions are quantities which are 
often measured in experiments. In [l| we investigated 
second order correlation functions and their remarkable 
differences between a canonical and a grand canonical 
description. We also showed the good agreement of 
our values with corrected grand canonical results from 
fl9j . Here we want to go beyond this well known re- 
sults and show correlation functions for the potential 




FIG. 5. First and second order spatial correlations Gi(a;,0) 
(top) and G2(x, 0) (bottom) obtained with our stochastic field 
equation for a Bose gas trapped in the potential V(x,y, z) = 
x 4 + y 4 + z 4 for different temperatures (0.5T C : crosses, 0.8T C : 
stars) . 



SFE harmonic 
■ analytical harmonic 




T/T c 



FIG. 6. Simulation (plus signs) of the variance of the ground 
state occupation for a 3D ideal Bose gas in an isotropic har- 
monic potential compared to results from [18( (full line). 



V(x,y,z) = x 4 + y 4 + z 4 . In Fig. [5] we present calcu- 
lations of first Gi(x, 0) = (i/>t (x)-0(O)) and second order 
G2(x,0) = (ip^ (x)^ (0)^(^)^(0)) correlation functions 
as they depend on the x-coordinate for a Bose gas of 200 
particles for temperatures below T c (plus signs, crosses). 

The differences of the grand canonical and the canon- 
ical ensemble become very obvious when considering the 
fluctuations of the ground state occupation as a function 
of temperature. As already mentioned in the introduc- 
tion, in a grand canonical description the variance for 
temperature near zero would be of the order of the par- 
ticle number TV. In Fig. [S] we display the ground state 
number fluctuations for a Bose gas of 200 particles in 
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an isotropic harmonic potential in the canonical ensem- 
ble that clearly tend to zero as temperature lowers. Wc 
compare our calculation with the analytical canonical re- 
sult of [lH which is based on the known eigenenergies for 
the harmonic potential. 

So far, for the determination of the graphs in Figs. [3K1 
is was sufficient to simulate with a single norm as elabo- 
rated upon in Sect. IIIII Now, we find that the number 
fluctuations are very sensitive with respect to slight er- 
rors of first (No) or second order (Nq) expectation values. 
Therefore, for this application, in order to reach the pre- 
cision shown in Fig. [5J it is necessary to simulate with 
a distribution of norms. We use the simplified Poisson- 
type distribution from (|T5j) . which explains why near the 
critical temperature deviations from the exact result oc- 
cur. 

One can summarize our simulations by pointing out 
that with our stochastic field equation (O in combina- 
tion with the Wigner-Weyl approximation for the noise 
generation, one gets good numerical results for many dif- 
ferent temperature dependent quantities in arbitrary po- 
tentials, also embracing the fixed particle number in these 
systems. 



VI. CONCLUSION AND OUTLOOK 

In this paper we have discussed a novel stochastic field 
equation for a Bose gas with a finite and fixed particle 
number: our theory is based on the canonical ensemble 
referring to the conditions found in actual experiments 
with atomic gases in traps. The equation is exact for 
ideal gases. A generalization including atomic interac- 
tions appears possible on the basis of the inclusion of a 
mean-field interaction term in the potential. 

Our approach to a numerical solution of the equation 
is explained in great detail. We focus on the implementa- 
tion in position space, in which it is not necessary to know 
the single-particle eigenenergies of the system. By chang- 
ing the specification of the external potential in the nu- 
merical code, we can easily determine equilibrium proper- 
ties of the gas for many trap geometries. We can calculate 
correlation functions of arbitrary order and in principle 
obtain information about the full thermal canonical state. 
We show that it is possible to achieve results with satisfy- 
ing precision for many different quantities of the ideal gas 
like first and second order correlation functions, ground 
state occupations and the variance of the ground state 
occupation. 
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Appendix A: Calculation of the norm distribution 

P(N) 



In this appendix we derive expression (|I3|) for the dis- 
tribution of the norm, 



P(A0 



J dfi{z}S(Af - \zk\ 2 )W N -i({z}) 
k 

/ dti{z}W N ({z}) 



; (Al) 



with 



N-l 



W N -l{z} 



(N-l) 



1M U> |S 



(A2) 

which we use in Section [III1 Suppose the number of con- 
sidered cigenstatcs equals M, which is arbitrary. We per- 
form the transformation of variables: \zi\ 2 — > X{ and find 

oo oo 

P(A0 = ly d M -Vx J dx 6(xo-(N--J2^)) 



1 



(N-l) 



(N-l) 



C(N — I)! 



jM 



_j - J2 (e^.-e^o), 



(A3) 



M-l 

0< Yl Xi<N 



:=F A 



with C = J dfi{z}WN({z}). First, we calculate the inte- 
gral F M - 



AT 

F M = / dxi 



dxz... 



dx 



M-l M — 1 

0< J2 Xi<N-X! 0<X M <N- J2 x i 

2=2 2 = 1 

M 

- E (e ^ -e f3e °)x i 
xe i=1 



(A4) 



Now we substitute xi = xi, X2 = xi + x\, xs = xs + x 2 
etc. and rewrite the integral as 

M AT 

F M = J dx ie -^ ei - e " 2 ^ J dx 2 e-^' 2 - e ' 3 ' 3 ^... 



o 



dx M e-^ M ~^°^. 



(A5) 



We are grateful for inspiring discussions with Markus 
Oberthaler and Thimo Grotz. Sigmund Heller ac- 
knowledges support from the International Max Planck 
Research School for Dynamical Processes in Atoms, 
Molecules and Solids, Dresden. 



By considering a derivative with respect to M of equa- 
tion (|A5|) it can be seen that Fm satisfies the differential 
equation 

ofm^o = e - (e ^. e ^o WFM _ im (A6) 
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We succeeded to find a solution with the ansatz bit of analysis we find 

M M 



k -° 1=0 
l^k 

M 

Fm(M) = TT — Finally, we plug expression (|A8|) into equation (|A3 

fe=l 
M 

+ E 



xx e /3 £fc _ e /3e ^ c = J-d/i{«}Wiv({«}) = J dAf J dfl{z}S(Af 

|zfe| 2 )M / Ar({z}) to write the norm distribution in the 



M M 

9k „-(e ge fc-e^°W fe 



-e 

fc— 1 



form 



(A7) AA^-DE^e- 



P(Af) 



k 



(N-l)\ J2c k e-^ N +^ 

k 



M-l M-l 

with < = n ^d .gf = _ e ,% +e ,. M . In this "° 

fc=l 

way we can obtain all coefficients <?j^ iteratively. With a (fT5|) of Section IIIII 



(A9) 

with c/c = ■^773^^-, as it is written in equation 
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